## author:    A. D?r, Robert A. Huber, Gemma Mateo, Gabriele Spilker
## contact:   robert.huber@sbg.ac.at
## file name: nti_regression.R
## Context:   Project on NTI in PTAs
## started:   2020-11-19
## Summary:   Provides main analyses of the paper

#https://cran.r-project.org/web/packages/clubSandwich/vignettes/panel-data-CRVE.html

# Function to extract quantities of interest for texreg 
source("./rcode/prepTex.R")
source("./rcode/prepTexglm.R")

m1r <- lm(rating ~ ser + ipr + env + lab, df_exp)

m1c <- glm(choice ~ ser + ipr + env + lab, df_exp, family = "binomial")

m1r5 <- glm(rating > 5 ~ ser + ipr + env + lab, df_exp, family = "binomial")

m2r <- lm(rating ~ ser + ipr + env + lab + ig, df_exp)

m2c <- glm(choice ~ ser + ipr + env + lab + ig, df_exp, family = "binomial")

m2r5 <- glm(I(rating > 5) ~ ser + ipr + env + lab + ig, df_exp, family = "binomial")

m3r <- lm(rating ~ (ser + ipr + env + lab) * ig, df_exp)

m3c <- glm(choice ~ (ser + ipr + env + lab) * ig, df_exp, family = "binomial")

m3r5 <- glm(I(rating > 5) ~ (ser + ipr + env + lab) * ig, df_exp, family = "binomial")

m4r <- lm(rating ~ (ser + ipr + env + lab) * ig_export, df_exp)

m4c <- glm(choice ~ (ser + ipr + env + lab) * ig_export, df_exp, family = "binomial")

m4r5 <- glm(I(rating > 5) ~ (ser + ipr + env + lab) * ig_export, df_exp, family = "binomial")

m5r <- lm(rating ~ ser * ig_ser + ipr + env + lab, df_exp)

m5c <- glm(choice ~ ser * ig_ser + ipr + env + lab, df_exp, family = "binomial")

m5r5 <- glm(I(rating > 5) ~ ser * ig_ser + ipr + env + lab, df_exp, family = "binomial")

m6r <- lm(rating ~ ser + ipr * ig_ipr + env + lab, df_exp)

m6c <- glm(choice ~ ser + ipr * ig_ipr + env + lab, df_exp, family = "binomial")

m6r5 <- glm(I(rating > 5) ~ ser + ipr * ig_ipr + env + lab, df_exp, family = "binomial")

m7r <- lm(rating ~ ser + ipr + env * ig_env + lab, df_exp)

m7c <- glm(choice ~ ser + ipr + env * ig_env + lab, df_exp, family = "binomial")

m7r5 <- glm(I(rating > 5) ~ ser + ipr + env * ig_env + lab, df_exp, family = "binomial")

m8r <- lm(rating ~ ser + ipr + env + lab * ig_lab, df_exp)

m8c <- glm(choice ~ ser + ipr + env + lab * ig_lab, df_exp, family = "binomial")

m8r5 <- glm(I(rating > 5) ~ ser + ipr + env + lab * ig_lab, df_exp, family = "binomial")

# Regression Tables ####

texreg::screenreg(list(prepTex(m1r), prepTex(m2r), prepTex(m3r)))

texreg::texreg(list(prepTex(m1r), prepTex(m2r), prepTex(m3r)),
               file = "./output/TableA4.tex",
               single.row = T,
               caption = "Design features, actors, and PTA rating",
               caption.above = T,
               label = "tab:regression",
               float.pos = "htb",
               custom.model.names = c("Model 1R", "Model 2R", "Model 3R"),
               custom.note = "\\parbox{0.8\\linewidth}{%stars. Estimates from an OLS regression. Standard errors in parentheses. Dependent variable (PTA rating) ranges from 1 to 10.}",
               custom.coef.names = c(NA,
                                     "Services", "IPR", "Environment", "Labour",
                                     "Export-oriented business",
                                     "Import-competing business",
                                     "Citizen Groups", "Labour unions",
                                     "Services x Export-oriented business",
                                     "Services x Import-competing business",
                                     "Services x Citizen groups",
                                     "Services x Labour unions",
                                     "IPR x Export-oriented business",
                                     "IPR x Import-competing business",
                                     "IPR x Citizen groups",
                                     "IPR x Labour unions",
                                     "Environ. x Export-oriented business",
                                     "Environ. x Import-competing business",
                                     "Environ. x Citizen groups",
                                     "Environ. x Labour unions",
                                     "Labour x Export-oriented business",
                                     "Labour x Import-competing business",
                                     "Labour x Citizen groups",
                                     "Labour x Labour unions"))

texreg::screenreg(list(prepTex(m4r), prepTex(m5r), prepTex(m6r), prepTex(m7r), prepTex(m8r)))

texreg::texreg(list(prepTex(m4r), prepTex(m5r), prepTex(m6r)),
               file = "./output/TableA5.tex",
               single.row = T,
               caption = "Additional evidence: specific actors, service and IPR provisions, and PTA rating",
               caption.above = T,
               label = "tab:additional_evidence1",
               float.pos = "htb",
               custom.model.names = c("Model 4R", "Model 5R", "Model 6R"),
               custom.note = "\\parbox{0.9\\linewidth}{%stars. Estimates from an OLS regression. Standard errors in parentheses. Dependent variable (PTA rating) ranges from 1 to 10.}",
               custom.coef.names = c(NA,
                                     "Services", "IPR", "Environment", "Labour",
                                     "Northern export-oriented business",
                                     "Service x Northern export-oriented business",
                                     "IPR x Northern export-oriented business",
                                     "Environ. x Northern export-oriented business",
                                     "Labour x Northern export-oriented business",
                                     "Service business",
                                     "Services x Service business",
                                     "Other citizen groups",
                                     "Knowledge intensive business",
                                     "Health policy citizen groups",
                                     "IPR x Other citizen groups",
                                     "IPR x Knowledge intensive business",
                                     "IPR x Health citizen groups"))

texreg::texreg(list(prepTex(m7r), prepTex(m8r)),
               file = "./output/TableA6.tex",
               single.row = T,
               caption = "Additional evidence: specific actors, environmental and labour provisions, and PTA rating",
               caption.above = T,
               label = "tab:additional_evidence2",
               float.pos = "htb",
               custom.model.names = c("Model 7R", "Model 8R"),
               custom.note = "\\parbox{.65\\linewidth}{%stars. Estimates from an OLS regression. Standard errors in parentheses. Dependent variable (PTA rating) ranges from 1 to 10.}",
               custom.coef.names = c(NA,
                                     "Services", "IPR", "Environment", 
                                     "Environmental citizen groups",
                                     "Labour",
                                     "Environ. x Environ. citizen groups",
                                     "Labour policy citizen groups",
                                     "Labour x Labour policy citizen groups"),
               reorder.coef = c(1:4, 6, 5, 7:9))

# ... Robustness checks ####

texreg::texreg(list(prepTexglm(m1r5), prepTexglm(m2r5), prepTexglm(m3r5)),
               file = "./output/TableA10.tex",
               single.row = T,
               caption = "Design features, actors, and PTA support (Rating above 5)",
               caption.above = T,
               label = "tab:regression_above5",
               float.pos = "htb",
               custom.model.names = c("Model 1S", "Model 2S", "Model 3S"),
               custom.note = "\\parbox{0.8\\linewidth}{%stars. Estimates from an logistic regression. Standard errors in parentheses. Dependent variable (PTA support) is a dummy.}",
               custom.coef.names = c(NA,
                                     "Services", "IPR", "Environment", "Labour",
                                     "Export-oriented business",
                                     "Import-competing business",
                                     "Citizen Groups", "Labour unions",
                                     "Services x Export-oriented business",
                                     "Services x Import-competing business",
                                     "Services x Citizen groups",
                                     "Services x Labour unions",
                                     "IPR x Export-oriented business",
                                     "IPR x Import-competing business",
                                     "IPR x Citizen groups",
                                     "IPR x Labour unions",
                                     "Environ. x Export-oriented business",
                                     "Environ. x Import-competing business",
                                     "Environ. x Citizen groups",
                                     "Environ. x Labour unions",
                                     "Labour x Export-oriented business",
                                     "Labour x Import-competing business",
                                     "Labour x Citizen groups",
                                     "Labour x Labour unions"))

texreg::texreg(list(prepTexglm(m4r5), prepTexglm(m5r5), prepTexglm(m6r5)),
               file = "./output/TableA11.tex",
               single.row = T,
               caption = "Additional evidence: specific actors, service and IPR provisions, and PTA support (Rating above 5)",
               caption.above = T,
               label = "tab:additional_evidence1_above5",
               float.pos = "htb",
               custom.model.names = c("Model 4S", "Model 5S", "Model 6S"),
               custom.note = "\\parbox{0.9\\linewidth}{%stars. Estimates from an logistic regression. Standard errors in parentheses. Dependent variable (PTA support) is a dummy.}",
               custom.coef.names = c(NA,
                                     "Services", "IPR", "Environment", "Labour",
                                     "Northern export-oriented business",
                                     "Service x Northern export-oriented business",
                                     "IPR x Northern export-oriented business",
                                     "Environ. x Northern export-oriented business",
                                     "Labour x Northern export-oriented business",
                                     "Service business",
                                     "Services x Service business",
                                     "Other citizen groups",
                                     "Knowledge intensive business",
                                     "Health policy citizen groups",
                                     "IPR x Other citizen groups",
                                     "IPR x Knowledge intensive business",
                                     "IPR x Health citizen groups"))

texreg::texreg(list(prepTexglm(m7r5), prepTexglm(m8r5)),
               file = "./output/TableA12.tex",
               single.row = T,
               caption = "Additional evidence: specific actors, environmental and labour provisions, and PTA support (Rating above 5)",
               caption.above = T,
               label = "tab:additional_evidence2_above5",
               float.pos = "htb",
               custom.model.names = c("Model 7S", "Model 8S"),
               custom.note = "\\parbox{.65\\linewidth}{%stars. Estimates from an logistic regression. Standard errors in parentheses. Dependent variable (PTA support) is a dummy.}",
               custom.coef.names = c(NA,
                                     "Services", "IPR", "Environment", 
                                     "Environmental citizen groups",
                                     "Labour",
                                     "Environ. x Environ. citizen groups",
                                     "Labour policy citizen groups",
                                     "Labour x Labour policy citizen groups"),
               reorder.coef = c(1:4, 6, 5, 7:9))


# ... Choice ####

texreg::texreg(list(prepTexglm(m1c), prepTexglm(m2c), prepTexglm(m3c)),
               file = "./output/TableA7.tex",
               single.row = T,
               caption = "Design features, actors, and PTA choice",
               caption.above = T,
               label = "tab:regression_choice",
               float.pos = "htb",
               custom.model.names = c("Model 1C", "Model 2C", "Model 3C"),
               custom.note = "\\parbox{0.8\\linewidth}{%stars. Estimates from an logistic regression. Standard errors in parentheses. Dependent variable (PTA choice) is a dummy.}",
               custom.coef.names = c(NA,
                                     "Services", "IPR", "Environment", "Labour",
                                     "Export-oriented business",
                                     "Import-competing business",
                                     "Citizen Groups", "Labour unions",
                                     "Services x Export-oriented business",
                                     "Services x Import-competing business",
                                     "Services x Citizen groups",
                                     "Services x Labour unions",
                                     "IPR x Export-oriented business",
                                     "IPR x Import-competing business",
                                     "IPR x Citizen groups",
                                     "IPR x Labour unions",
                                     "Environ. x Export-oriented business",
                                     "Environ. x Import-competing business",
                                     "Environ. x Citizen groups",
                                     "Environ. x Labour unions",
                                     "Labour x Export-oriented business",
                                     "Labour x Import-competing business",
                                     "Labour x Citizen groups",
                                     "Labour x Labour unions"))



texreg::texreg(list(prepTexglm(m4c), prepTexglm(m5c), prepTexglm(m6c)),
               file = "./output/TableA8.tex",
               single.row = T,
               caption = "Additional evidence: specific actors, service and IPR provisions, and PTA choice",
               caption.above = T,
               label = "tab:additional_evidence1_choice",
               float.pos = "htb",
               custom.model.names = c("Model 4C", "Model 5C", "Model 6C"),
               custom.note = "\\parbox{.9\\linewidth}{%stars. Estimates from an logistic regression. Standard errors in parentheses. Dependent variable (PTA choice) is a dummy.}",
               custom.coef.names = c(NA,
                                     "Services", "IPR", "Environment", "Labour",
                                     "Northern export-oriented business",
                                     "Service x Northern export-oriented business",
                                     "IPR x Northern export-oriented business",
                                     "Environ. x Northern export-oriented business",
                                     "Labour x Northern export-oriented business",
                                     "Service business",
                                     "Services x Service business",
                                     "Other citizen groups",
                                     "Knowledge intensive business",
                                     "Health policy citizen groups",
                                     "IPR x Other citizen groups",
                                     "IPR x Knowledge intensive business",
                                     "IPR x Health citizen groups"))

texreg::texreg(list(prepTexglm(m7c), prepTexglm(m8c)),
               file = "./output/TableA9.tex",
               single.row = T,
               caption = "Additional evidence: specific actors, environmental and labour provisions, and PTA choice",
               caption.above = T,
               label = "tab:additional_evidence2_choice",
               float.pos = "htb",
               custom.model.names = c("Model 7C", "Model 8C"),
               custom.note = "\\parbox{0.65\\linewidth}{%stars. Estimates from an logistic regression. Standard errors in parentheses. Dependent variable (PTA choice) is a dummy.}",
               custom.coef.names = c(NA,
                                     "Services", "IPR", "Environment", 
                                     "Environmental citizen groups",
                                     "Labour",
                                     "Environ. x Environ. citizen groups",
                                     "Labour policy citizen groups",
                                     "Labour x Labour policy citizen groups"),
               reorder.coef = c(1:4, 6, 5, 7:9))

# Figures ####

summary(margins::margins(m3r,
                         variables = c("ser", "ipr", "env", "lab"),
                         vcov = clubSandwich::vcovCR(m3r, type = "CR0", cluster = df_exp$respondent),
                         at = list(ig = c("Other business", "Export oriented business", "Import competing business","Citizen groups", "Labour unions")))) %>%
#c("Southern import-competing business", "Southern export-oriented business", "Northern import-competing business", "Northern export-oriented business", "Citizen groups", "Labour unions")))) %>%
  mutate(lower95 = lower,
         upper95 = upper,
         lower90 = AME - 1.645 * SE,
         upper90 = AME + 1.645 * SE,
         feature = factor(ifelse(grepl(pattern = "ser", factor), "Services",
                                 ifelse(grepl(pattern = "IPR", factor), "IPR",
                                        ifelse(grepl(pattern = "env", factor), "Environment",
                                               ifelse(grepl(pattern = "lab", factor), "Labour rights", NA)))),
                          levels = c("Services", "IPR", "Environment", "Labour rights")),
         ig = factor(ig,
                     levels = c("Other business", "Export oriented business", "Import competing business","Citizen groups", "Labour unions"),
                     labels = c("Other business", "Export-oriented business", "Import-competing business","Citizen groups", "Labour unions"))) %>%
  ggplot(., aes(forcats::fct_rev(ig), AME)) + 
  geom_hline(yintercept = 0, lty = "dashed") + 
  geom_linerange(aes(ymin = lower90, ymax = upper90), size = 1.5, position = position_dodge(width = 1)) +
  geom_pointrange(aes(ymin = lower95, ymax = upper95), position = position_dodge(width = 1)) +
  #  ggtitle("Service Provisions") + 
  facet_grid(feature~.) + 
  theme_minimal() +
  coord_flip() + 
  theme(panel.border = element_rect(color = "black", fill = NA, size = 1),
        legend.position = "none") +
  labs(x = "", y = "\nConditional Average Marginal Component\nEffects of Including Provisions")

ggsave("./output/Figure2.eps",
       device = "eps", 
       width = 15, height = 12, units = "cm")

summary(margins::margins(m4r,
                         variables = c("ser", "ipr", "env", "lab"),
                         vcov = clubSandwich::vcovCR(m4r, type = "CR0", cluster = df_exp$respondent),
                         at = list(ig_export = c("Southern export oriented business", "Northern export oriented business")))) %>%
  mutate(lower95 = lower,
         upper95 = upper,
         lower90 = AME - 1.645 * SE,
         upper90 = AME + 1.645 * SE,
         feature = factor(ifelse(grepl(pattern = "ser", factor), "Services",
                                 ifelse(grepl(pattern = "IPR", factor), "IPR",
                                        ifelse(grepl(pattern = "env", factor), "Environment",
                                               ifelse(grepl(pattern = "lab", factor), "Labour rights", NA)))),
                          levels = c("Services", "IPR", "Environment", "Labour rights"),
                          labels = c("Services", "IPR", "Environ.", "Labour")),
         ig_export = factor(ig_export,
                            levels = c("Southern export oriented business", "Northern export oriented business"),
                            labels = c("Southern export-oriented business", "Northern export-oriented business"))) %>%
  ggplot(., aes(ig_export, AME)) + 
  geom_hline(yintercept = 0, lty = "dashed") + 
  geom_linerange(aes(ymin = lower90, ymax = upper90), size = 1.5, position = position_dodge(width = 1)) +
  geom_pointrange(aes(ymin = lower95, ymax = upper95), position = position_dodge(width = 1)) +
  #  ggtitle("Service Provisions") + 
  facet_grid(feature~.) + 
  theme_minimal() +
  coord_flip() + 
  theme(panel.border = element_rect(color = "black", fill = NA, size = 1),
        legend.position = "none") +
  labs(x = "", y = "\nConditional Average Marginal Component\nEffects of Including Provisions")

ggsave("./output/Figure4.eps",
       device = "eps",
       width = 15, height = 8, units = "cm")


bind_rows(summary(margins::margins(m5r,
                                   variables = c("ser"),
                                   vcov = clubSandwich::vcovCR(m5r, type = "CR0", cluster = df_exp$respondent),
                                   at = list(ig_ser = c("Other business", "Service business")))) %>%
            dplyr::mutate(ig = ig_ser) %>%
            select(-ig_ser),
          summary(margins::margins(m6r,
                                   variables = c("ipr"),
                                   vcov = clubSandwich::vcovCR(m6r, type = "CR0", cluster = df_exp$respondent),
                                   at = list(ig_ipr = c("Other business", "Other citizen group", "Knowledge intensive business", "Health policy citizen group")))) %>%
            mutate(ig = ig_ipr) %>%
            select(-ig_ipr),
          summary(margins::margins(m7r,
                                   variables = c("env"),
                                   vcov = clubSandwich::vcovCR(m7r, type = "CR0", cluster = df_exp$respondent),
                                   at = list(ig_env = c("Other citizen group", "Environmental citizen group")))) %>%
            mutate(ig = ig_env) %>%
            select(-ig_env),
          summary(margins::margins(m8r,
                                   variables = c("lab"),
                                   vcov = clubSandwich::vcovCR(m8r, type = "CR0", cluster = df_exp$respondent),
                                   at = list(ig_lab = c("Other citizen group", "Labour policy citizen group")))) %>%
            mutate(ig = ig_lab) %>%
            select(-ig_lab)) %>%
  mutate(lower95 = lower,
         upper95 = upper,
         lower90 = AME - 1.645 * SE,
         upper90 = AME + 1.645 * SE,
         feature = factor(ifelse(grepl(pattern = "ser", factor), "Services",
                                 ifelse(grepl(pattern = "IPR", factor), "IPR",
                                        ifelse(grepl(pattern = "env", factor), "Environment",
                                               ifelse(grepl(pattern = "lab", factor), "Labour rights", NA)))),
                          levels = c("Services", "IPR", "Environment", "Labour rights"),
                          labels = c("Services", "IPR", "Environ.", "Labour")),
         ig = factor(ig,
                     levels = c("Service business",
                                "Knowledge intensive business",
                                "Other business",
                                "Health policy citizen group",
                                "Other citizen group",
                                "Labour policy citizen group",
                                "Environmental citizen group"),
                     labels = c("Tradeable service",
                                "Knowledge-intensive business",
                                "Rest of business",
                                "Health policy citizen group",
                                "Other citizen group",
                                "Labour policy citizen group",
                                "Environmental citizen group"))) %>%
  ggplot(., aes(forcats::fct_rev(ig), AME)) + 
  geom_hline(yintercept = 0, lty = "dashed") + 
  geom_linerange(aes(ymin = lower90, ymax = upper90), size = 1.5, position = position_dodge(width = 1)) +
  geom_pointrange(aes(ymin = lower95, ymax = upper95), position = position_dodge(width = 1)) +
  #  ggtitle("Service Provisions") + 
  facet_grid(feature~., scales = "free", space = "free") + 
  theme_minimal() +
  coord_flip() + 
  theme(panel.border = element_rect(color = "black", fill = NA, size = 1),
        legend.position = "none") +
  scale_y_continuous(breaks = seq(from = -2, to = 4, by = 1)) + 
  labs(x = "", y = "\nConditional Average Marginal Component\nEffects of Including Provisions")

ggsave("./output/Figure5.eps",
       device = "eps", 
       width = 15, height = 10, units = "cm")


# Scenarios ####
pred <- data.frame(ser = rep(c("Services:No", "Services:Yes"),2),
                   ipr = rep(c("IPR:No", "IPR:Yes"),2),
                   env = rep(c("Environmental:No", "Environmental:Yes"), each = 2),
                   lab = rep(c("Labour:No", "Labour:Yes"), each = 2))

pred_in <- rbind(pred %>%
                   mutate(ig = "Other business"),
                 pred %>%
                   mutate(ig = "Export oriented business"),
                 pred %>%
                   mutate(ig = "Import competing business"),
                 pred %>%
                   mutate(ig = "Citizen groups"),
                 pred %>%
                   mutate(ig = "Labour unions")) %>%
  mutate(ig = factor(ig,
                     levels = c("Other business", "Export oriented business",
                                "Import competing business",
                                "Citizen groups", "Labour unions")))

pred_out <- left_join(pred_in, prediction::prediction_summary(model = m3r,
                                                              at = list(ser = c("Services:No", "Services:Yes"),
                                                                        ipr = c("IPR:No", "IPR:Yes"),
                                                                        env = c("Environmental:No", "Environmental:Yes"),
                                                                        lab = c("Labour:No", "Labour:Yes"),
                                                                        ig = c("Other business", "Export oriented business",
                                                                               "Import competing business",
                                                                               "Citizen groups", "Labour unions")),
                                                              vcov = clubSandwich::vcovCR(m3r, type = "CR0", cluster = df_exp$respondent)),
                      by = c("ser" = "at(ser)", "ipr" = "at(ipr)",
                             "env" = "at(env)", "lab" = "at(lab)", "ig" = "at(ig)")) %>%
  mutate(lower95 = lower,
         upper95 = upper,
         lower90 = Prediction - 1.645 * SE,
         upper90 = Prediction + 1.645 * SE,
         agreement = factor(rep(c("Empty", "Trade", "Non-Trade", "Full"), 5),
                            levels = c("Empty", "Trade", "Non-Trade", "Full"),
                            labels = c("Bare-bones", "Trade", "Non-Trade", "Comprehensive")),
         ig = factor(ig,
                     levels = c("Other business", "Export oriented business",
                                "Import competing business",
                                "Citizen groups", "Labour unions"),
                     labels = c("Other\nbusiness", "Export\noriented\nbusiness",
                                "Import\ncompeting\nbusiness",
                                "Citizen\ngroups", "Labour\nunions")))

ggplot(pred_out, aes(x = forcats::fct_rev(agreement),
                 y = Prediction)) + 
  geom_hline(yintercept = 5.5, lty = "dotted") + 
  geom_pointrange(aes(ymin = lower95,
                      ymax = upper95)) + 
  geom_linerange(aes(ymin = lower90,
                     ymax = upper90),
                 size = 1.5) +
  facet_grid(ig~.) + 
  labs(x = "Scenario", y = "Predicted Rating") + 
  theme_minimal() +
  coord_flip() + 
  theme(panel.border = element_rect(color = "black", fill = NA, size = 1),
        legend.position = "none") +
  scale_y_continuous(breaks = c(1:7)) + 
  NULL

ggsave("./output/Figure3.eps",
       device = "eps",
       width = 15, height = 10, units = "cm")

# Leeper ####
m1r_leeper <- rating ~ ser + ipr + env + lab

cregg::cj(df_exp, m1r_leeper, id = ~ respondent, estimate = "mm",  h0 = 5.5, by = ~ig) %>%
  mutate(lower95 = lower,
         upper95 = upper,
         lower90 = estimate - 1.645 * std.error,
         upper90 = estimate + 1.645 * std.error,
         factor = factor(x = feature,
                         levels = c("ser", "ipr", "env", "lab"),
                         labels = c("Services", "IPR", "Environment", "Labour rights")),
         dim = factor(x = ifelse(grepl("No", x = level), "Excluded",
                                 ifelse(grepl("Yes", x = level), "Included", NA)),
                      levels = c("Excluded", "Included"))) %>%
  ggplot(., aes(x = fct_rev(ig), y = estimate, colour = dim, shape = dim)) + 
  geom_hline(yintercept = 5.5, lty = "dotted") + 
  geom_pointrange(aes(ymin = lower95,
                      ymax = upper95), 
                  position = position_dodge(width = 1)) + 
  geom_linerange(aes(ymin = lower90,
                     ymax = upper90),
                 size = 1.5,
                 position = position_dodge(width = 1)) +
  facet_grid(factor~., scales = "free_y") + 
  coord_flip() +
  theme_minimal() +
  labs(x = "", y = "Estimate") + 
  scale_colour_discrete("Group type") + 
  scale_shape_discrete("Group type") +
  theme(panel.border = element_rect(color = "black", fill = NA, size = 1),
        legend.position = "bottom") +
  NULL

ggsave("./output/FigureA2.pdf",
       width = 15,
       height = 13,
       units = "cm")

cregg::cj(df_exp, m1r_leeper, id = ~ respondent, estimate = "mm",  h0 = 5.5, by = ~ig_export) %>%
  mutate(lower95 = lower,
         upper95 = upper,
         lower90 = estimate - 1.645 * std.error,
         upper90 = estimate + 1.645 * std.error,
         factor = factor(x = feature,
                         levels = c("ser", "ipr", "env", "lab"),
                         labels = c("Services", "IPR", "Environ.", "Labour")),
         dim = factor(x = ifelse(grepl("No", x = level), "Excluded",
                                 ifelse(grepl("Yes", x = level), "Included", NA)),
                      levels = c("Excluded", "Included"))) %>%
  ggplot(., aes(x = ig_export, y = estimate, colour = dim, shape = dim)) + 
  geom_hline(yintercept = 5.5, lty = "dotted") + 
  geom_pointrange(aes(ymin = lower95,
                      ymax = upper95), 
                  position = position_dodge(width = 1)) + 
  geom_linerange(aes(ymin = lower90,
                     ymax = upper90),
                 size = 1.5,
                 position = position_dodge(width = 1)) +
  facet_grid(factor~., scales = "free_y") + 
  coord_flip() +
  theme_minimal() +
  labs(x = "", y = "Estimate") + 
  scale_colour_discrete("Provision") + 
  scale_shape_discrete("Provision") +
  theme(panel.border = element_rect(color = "black", fill = NA, size = 1),
        legend.position = "bottom") +
  NULL

ggsave("./output/FigureA3.pdf",
       width = 15,
       height = 10,
       units = "cm")


bind_rows(cregg::cj(df_exp, m1r_leeper, id = ~ respondent, estimate = "mm",  h0 = 5.5, by = ~ig_ser) %>%
            filter(feature == "ser") %>%
            select(ig = ig_ser, everything()),
          cregg::cj(df_exp, m1r_leeper, id = ~ respondent, estimate = "mm",  h0 = 5.5, by = ~ig_ipr) %>%
            filter(feature == "ipr") %>%
            select(ig = ig_ipr, everything()),
          cregg::cj(df_exp, m1r_leeper, id = ~ respondent, estimate = "mm",  h0 = 5.5, by = ~ig_env) %>%
            filter(feature == "env") %>%
            select(ig = ig_env, everything()),
          cregg::cj(df_exp, m1r_leeper, id = ~ respondent, estimate = "mm",  h0 = 5.5, by = ~ig_lab) %>%
            filter(feature == "lab") %>%
            select(ig = ig_lab, everything())) %>%
            mutate(lower95 = lower,
                   upper95 = upper,
                   lower90 = estimate - 1.645 * std.error,
                   upper90 = estimate + 1.645 * std.error,
                   factor = factor(x = feature,
                                   levels = c("ser", "ipr", "env", "lab"),
                           labels = c("Services", "IPR", "Environ.", "Labour")),
                   dim = factor(x = ifelse(grepl("No", x = level), "Excluded",
                                           ifelse(grepl("Yes", x = level), "Included", NA)),
                                levels = c("Excluded", "Included")),
                   ig = factor(ig,
                               levels = c("Service business",
                                          "Knowledge intensive business",
                                          "Other business",
                                          "Health policy citizen group",
                                          "Other citizen group",
                                          "Labour policy citizen group",
                                          "Environmental citizen group"),
                               labels = c("Tradeable service",
                                          "Knowledge-intensive business",
                                          "Rest of business",
                                          "Health policy citizen group",
                                          "Other citizen group",
                                          "Labour policy citizen group",
                                          "Environmental citizen group"))) %>%
  ggplot(., aes(forcats::fct_rev(ig), estimate, colour = dim, shape = dim)) + 
  geom_hline(yintercept = 5.5, lty = "dashed") + 
  geom_linerange(aes(ymin = lower90, ymax = upper90), size = 1.5, position = position_dodge(width = 1)) +
  geom_pointrange(aes(ymin = lower95, ymax = upper95), position = position_dodge(width = 1)) +
  #  ggtitle("Service Provisions") + 
  facet_grid(factor~., scales = "free", space = "free") + 
  theme_minimal() +
  coord_flip() + 
  theme(panel.border = element_rect(color = "black", fill = NA, size = 1),
        legend.position = "bottom") +
  scale_colour_discrete("Provision") + 
  scale_shape_discrete("Provision") +
  #scale_y_continuous(breaks = seq(from = -2, to = 4, by = 1)) + 
  labs(x = "", y = "\nConditional Average Marginal Component\nEffects of Including Provisions")

ggsave("./output/FigureA4.pdf", width = 15, height = 10, units = "cm")

